#coding=utf-8

import numpy as np
from numpy import linalg
from scipy.stats import f

#多元线性回归
class myMLR:
    def __init__(self,X,Y):
        self.X = X
        self.Y = Y
        self.A=None
    def fit(self):
        Xt= self.X.T
        XtX=np.dot(Xt,self.X)
        XtX_inv=np.linalg.inv(XtX)
        temp=np.dot(XtX_inv,Xt)
        self.A=np.dot(temp,self.Y)
    def getCoef(self):
        return self.A
    def predict(self,X):
        Y = np.dot(X,self.A)
        return Y
    #F检验
    def Ftest(self,alpha):
        n=len(self.X)
        k=self.X.shape[-1]-1
        f_arfa=f.isf(alpha,k,n-k-1)
        Yaver=self.Y.sum()/n
        Yhat=self.predict(self.X)
        U=((Yhat-Yaver)**2).sum()
        Qe=((self.Y-Yhat)**2).sum()
        F=(U/k)/(Qe/(n-k-1))
        answer=[F,f_arfa,F>f_arfa]
        return answer
#主成分分析
class myPCA:
    def __init__(self,X):
        self.X=X
        self.__lambda=None
    def SVD_decompose(self):
        B=np.linalg.svd(self.X,full_matrices=False)
        
        self.__lambda= lamda = B[1]
        sumtotual = np.sum(lamda)
        self.P=B[2].T
        self.T=B[0]*B[1]
        compare=[]
        sigmalist=[]
        sumv=0.0
        for i in range(len(lamda)-1):
            compare.append(lamda[i]/lamda[i+1])
            sumv+=lamda[i]/sumtotual
            sigmalist.append(sumv)
        # compare = [lamda[i]/lamda[i+1] for i in range(len(lamda)-1)]
        return np.array(compare),np.array(sigmalist)
    def PCAdecompose(self,k):
        T= self.T[:,:k]
        P=self.P[:,:k]
        return T,P
    def CallAll(self):
        compare=self.SVD_decompose()
        print("相邻特征值比值：",compare)
        k = int(input("请输入主成分数>>"))
        return self.PCAdecompose(k)
    def Auto_decompose(self,threshold=0.98):
        B=np.linalg.svd(self.X,full_matrices=False)
        self.__lambda= lamda = B[1]
        sumtotual = np.sum(lamda)
        self.P=B[2].T
        self.T=B[0]*B[1]
        sumv=0.0
        count = 0
        for i in range(len(lamda)-1):
            count+=1
            sumv+=lamda[i]/sumtotual
            if sumv>=threshold:
                break
        return self.PCAdecompose(count)
#主成分回归
class myPCR:
    def __init__(self,X,Y) -> None:
        self.X = X
        self.Y = Y
        self.pca=None
        self.T=None
        self.P=None
        self.mlr=None
        self.A=None
    def confirmPCs(self):
        self.pca=myPCA(self.X)
        compares = self.pca.SVD_decompose()
        return compares
    def fit(self,k):
        T,P=self.pca.PCAdecompose(k)
        self.P=P
        oneCol=np.ones(T.shape[0])
        Tnew=np.c_[oneCol,T]
        self.mlr=myMLR(Tnew,self.Y)
        self.mlr.fit()
        self.A=self.mlr.getCoef()
    def Auto_confirmfit(self,threshold=0.95):
        self.pca=myPCA(self.X)
        T,P=self.pca.Auto_decompose(threshold)
        self.P=P
        oneCol=np.ones(T.shape[0])
        Tnew=np.c_[oneCol,T]
        self.mlr=myMLR(Tnew,self.Y)
        self.mlr.fit()
        self.A=self.mlr.getCoef()
    def predict(self,Xnew):
        T = np.dot(Xnew,self.P)
        # print(T)
        # oneCol = np.ones(T.shape[0])
        T = np.r_[1,T]
        # T=np.c_[oneCol,T]
        ans=self.mlr.predict(T)
        return ans

def principal_comps(dataset, threshold = 0.98):
    # 返回满足要求的特征向量
    ret = []
    data = []

    # 标准化
    for (index, line) in enumerate(dataset):
        dataset[index] -= np.mean(line)
        # np.std(line, ddof = 1)即样本标准差(分母为n - 1)
        dataset[index] /= np.std(line, ddof = 1)
    # 求协方差矩阵
    Cov = np.cov(dataset)
    # 求特征值和特征向量
    eigs, vectors = linalg.eig(Cov)
    # 第i个特征向量是第i列，为了便于观察将其转置一下
    for i in range(len(eigs)):
        data.append((eigs[i], vectors[:, i].T))
    # 按照特征值从大到小排序
    data.sort(key = lambda x: x[0], reverse = True)

    sum = 0
    yi=1
    for comp in data:
        sum += comp[0] / np.sum(eigs)
        ret.append(
            tuple(map(
                # 保留5位小数
                lambda x: np.round(x, 5),
                # 特征向量、方差贡献率、累计方差贡献率
                (comp[1], comp[0] / np.sum(eigs), sum)
            ))
        )
        mylist = ret[-1][0]
        # print(to_latexf(mylist,yi)+"\\\\")
        yi+=1
        print('特征值:', comp[0], '特征向量:', ret[-1][0], '方差贡献率:', ret[-1][1], '累计方差贡献率:', ret[-1][2])
        if sum > threshold:
            return ret      

    return ret   
